
source("3.2.1 EtaEst.R")
source("3.2.2 GetRatio.R")
source("3.2.3 GetMaxKGen.R")
source("3.2.4 GetK.R")
source("3.2.5 PMaxEstRough.R")
source("3.2.6 nLL.R")
source("3.2.7 MarginalRatio.R")

source("3.5.1 DR.R")
source("3.5.2 PropScoreEst.R")
source("3.5.3 DR_MyFunc.R")
source("3.5.4 DR_EE_Value.R")

DREst = function(data, cnt, K, seed, max.step = 500, thres = 1e-4, 
                 mle, outcome.type = "MLE", dr.warm = NULL){
  # Output: theta, time, contrast
  # nll is just filled with 0---DR has nothing to do with nll
  # max.step = 5;thres = 1e-4; cnt <- count; outcome.type = "MLE" #DEBUG
  
  # Model for E[L_k+1|A_k, L_k, X]
  eta.coef        = EtaEst(data, cnt, K)
  X.names         = paste("X",1:pL0,sep="")
  names(eta.coef) = c("gamma00","gamma01","gamma10","gamma11","gamma.",X.names)
  
  # eta.coef        = matrix(eta.coef[c(1,6:(5+pL0), 2, 6:(5+pL0), 3, 6:(5+pL0),
  #                                     4, 6:(5+pL0), 5:(5+pL0))],,5)
  # colnames(eta.coef) = c("00","01","10","11",".")
  
  # Model for Propensity Score
  b.coef  = B.Est(data, cnt, K)
  
  # Model for the outcome
  p.y.func <- NA
  theta.coef = mle$theta.coef[c(1:2, 4:7)]
  # phi.coef = runif(pL0+3,  -epsilon, epsilon) # check dimension
  phi.coef = mle$phi.coef
  gop.coef = mle$gop.coef[c(1, 4:7)]
  source("3.5.3 DR_MyFunc.R")
  if (outcome.type=="MLE"){
    p.y.func <- get.y.prob.from.MLE(data, theta.coef, eta.coef,gop.coef,phi.coef,
                                    cnt, rep=500, seed)
  }else if(outcome.type=="logit"){
    p.y.func = y.logit.est(data, cnt, K)
  }else{
    stop("Outcome.type must be either MLE or logit")
  }
  
    
  set.seed(seed)
  ptm = proc.time()
  point.est <- dr.estimate.onestep(data, theta.coef, eta.coef,gop.coef, 
                                   p.y.func, b.coef, 
                                   cnt, rep=500, seed,
                                   max.step = max.step, thres = thres,
                                   dr.warm=dr.warm)
  time = (proc.time() - ptm)[3] 
  
  pars <- list()
  pars$theta.coef = point.est
  pars$gop.coef = gop.coef
  pars$phi.coef = mle$phi.coef
  eta.coef        = matrix(eta.coef[c(1,6:(5+pL0), 2, 6:(5+pL0), 3, 6:(5+pL0),
                                      4, 6:(5+pL0), 5:(5+pL0))],,5)
  colnames(eta.coef) = c("00","01","10","11",".")
  
  
  
  pars$contrast = CalMarRatio(pars, data, cnt, eta.coef, rep = 500,seed=seed)
  
  pars$theta.coef = c(pars$theta.coef[1:2],NA,pars$theta.coef[3:(pL0+2)])
  pars$gop.coef   = c(pars$gop.coef[1],NA,NA,pars$gop.coef[2:(pL0+1)])
  pars$eta.coef   = t(matrix(c(eta.coef[1,c(1:2,5)],eta.coef[2:(pL0+1),1],
                               eta.coef[1,3:4], rep(NA,pL0+1)),,2))
  pars$time = time 
  
  pars$n01  = NA
  pars$nll = DR.EE.Val(data, cnt, K, seed, max.step = max.step, thres =thres, 
                       mle,  outcome.type = outcome.type,dr.est=pars)
  
  return(pars)
  
}
